#ifndef _MS_DISTRIBUTION_H_

#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <fstream>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_statistics.h>

#include "fe_vecmat.h"
#include "ms_misc.h"

using std::ofstream;
/*
	Class structure to access the different
	random distributions provided by the GNU GSL
	
	And some statistical tools (histograms, ...)
*/

/* ------------------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------------------ */
/* Abstract class  -  serves as a container for all other members */
class RandomDistribution; 	/* abstract base class */
class UniformDistribution;	/*  */
class GaussianDistribution;	/*  */
class LogNormalDistribution; 	/*  */
class ExponentialDistribution; 	/*  */
class LaPlaceDistribution; 	/*  */
class CauchyDistribution;	/*  */
class RayleighDistribution;	/*  */
class LandauDistribution;	/*  */
class LevyAlphaStableDistribution;	/*  */
class LevySkewAlphaStableDistribution;	/*  */
class GammaDistribution;	/*  */
class ChiSquaredDistribution;	/*  */
class FDistribution;    	/*  */
class TDistribution;  /* Student's t Distribution */
class BetaDistribution; 	/*  */
class LogisticDistribution;	/*  */
class WeibullDistribution;	/*  */
class Gumbel1Distribution;	/* Typ 1 */
class Gumbel2Distribution;	/* Typ 2 */
// class DirichletDistribution;
class PoissonDistribution;	/*  */
class BinomialDistribution; 	/* Binomial? unsigned int n? TODO! */
class PascalDistribution;	/*  */
class LogarithmicDistribution;	/*  */


/* ------------------------------------------------------------------------------------ */
// class MisesFisherDistribution;
/* ------------------------------------------------------------------------------------ */
class Histogram;
/* ------------------------------------------------------------------------------------ */

/* ------------------------------------------------------------------------------------ */
class RandomDistribution
{
	protected:
		const gsl_rng_type *T;
		gsl_rng *r;
		// nachschauen wie man private nach unten vererben kann. friend? protected?
	public:
		RandomDistribution();
		void setSeed(unsigned long int s);
		virtual double	GetRandomValue() const = 0;	/* Return one random value for the current distribtion */
		virtual void	GetRandomVector(double * x, const int N ) const = 0;	/* Return N random values for the current distribtion */
		inline  void	GetRandomVector( Vector & x ) const
				{ GetRandomVector( x.X, x.n ); }	/* Return N random values for the current distribtion */
		double	RandVal() const;					/* return a uniformly distributed scalar value on [0, 1] */
		void		RandVal( double * x, const int N );		/* return a vector of uniformly distributed scalar value on [0, 1] */
		inline void		RandVal( Vector & x ) { RandVal( x.X, x.n ); }		/* return a vector of uniformly distributed scalar value on [0, 1] */
		virtual ~RandomDistribution();
};
/* ------------------------------------------------------------------------------------ */


class UniformDistribution : public RandomDistribution
{
	protected:
		double a,b;

	public:
		UniformDistribution( const double aa = 0., const double bb = 1. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
		inline double	GetParameters( double &aa, double &bb) const { aa = a; bb = b; return a; return b;}
};


class GaussianDistribution : public RandomDistribution
{
	protected:
		double sigma;

	public:
		
		GaussianDistribution( const double sigma = 1. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */

};


class LogNormalDistribution : public RandomDistribution
{
	protected:
		double zeta, sigma;

	public:
		
		LogNormalDistribution( const double zeta = 1., const double sigma = 1. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */

};


class ExponentialDistribution : public RandomDistribution
{
	protected:
		double mu;

	public:
		
		ExponentialDistribution( const double mu = 1.);
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */

};

class LaPlaceDistribution : public RandomDistribution
{
	protected:
		double a;

	public:
		
		LaPlaceDistribution( const double a = 1.);
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */

};


class CauchyDistribution  : public RandomDistribution
{
	protected:
		double a;

	public:
		
		CauchyDistribution( const double a = 0.);
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */

};


class RayleighDistribution : public RandomDistribution
{
	protected:
		double sigma;

	public:
		
		RayleighDistribution( const double sigma = 3. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class LandauDistribution : public RandomDistribution
{

	public:
		LandauDistribution();
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class LevyAlphaStableDistribution: public RandomDistribution
{
	protected:
		double c, alpha;

	public:
		
		LevyAlphaStableDistribution( const double c = 0., const double alpha = 1. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class LevySkewAlphaStableDistribution : public RandomDistribution
{
	protected:
		double c , alpha , beta;

	public:
		LevySkewAlphaStableDistribution( const double c = 0. , const double alpha = 1., const double beta = 2. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class GammaDistribution : public RandomDistribution
{
	protected:
		double a, b;

	public:
		
		GammaDistribution( const double a = 0., const double b = 1. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class ChiSquaredDistribution : public RandomDistribution
{
	protected:
		double nu;

	public:
		
		ChiSquaredDistribution( const double nu = 0. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class FDistribution : public RandomDistribution
{
	protected:
		double nu1, nu2;

	public:
		
		FDistribution( const double nu1 = 0., const double nu2 = 1. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class TDistribution : public RandomDistribution
{
	protected:
		double nu;

	public:
		
		TDistribution( const double nu = 0. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class BetaDistribution : public RandomDistribution
{
	protected:
		double a, b;

	public:
		
		BetaDistribution( const double a = 0., const double b  = 1. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class LogisticDistribution : public RandomDistribution
{
	protected:
		double a;

	public:
		
		LogisticDistribution( const double a = 0. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class WeibullDistribution : public RandomDistribution
{
	protected:
		double a, b;

	public:
		
		WeibullDistribution( const double a = 0., const double b = 1. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class Gumbel1Distribution : public RandomDistribution
{
	protected:
		double a , b;

	public:
		
		Gumbel1Distribution( const double a = 0., const double b = 1. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class Gumbel2Distribution : public RandomDistribution
{
	protected:
		double a , b; 

	public:
		
		Gumbel2Distribution( const double a = 0., const double b = 1. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};


class PoissonDistribution : public RandomDistribution
{
	protected:
		double mu;

	public:
		
		PoissonDistribution( const double mu = 0. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class BinomialDistribution : public RandomDistribution
{
	protected:
		double p;
		unsigned int n; 

	public:
		
		BinomialDistribution( const double p = 0., const unsigned int n = 1);
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
} ;

class PascalDistribution : public RandomDistribution
{
	protected:
		double p;
		unsigned int n;

	public:
		
		PascalDistribution( const double p = 0., const unsigned int n = 1);
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
};

class LogarithmicDistribution : public RandomDistribution
{
	protected:
		double p;

	public:
		
		LogarithmicDistribution( const double p = 0. );
		double	GetRandomValue() const ;
		void	GetRandomVector(double * x, const int N ) const ;	/* Return N random values for the current distribtion */
} ;


/* ------------------------------------------------------------------------------------ */
class MisesFisherDistribution
{
	protected:
	
	

	public:
		MisesFisherDistribution();
	
	
	
};




/* ------------------------------------------------------------------------------------ */


class Histogram {
	private:
		void	Recalculate();			/* Recalculate the histogram, e.g., if N changes or the histogram is initialized */
		gsl_histogram * gsl_h;
	protected:
		double	v_min, v_max, h_bin, SummeWW;		/* min/max value (interval boundaries), interval width */
		int	n_bin;				/* number of bins/intervals */
		Vector	hist, dat, W;			/* histogram data; data vector, Weighted Vector */
		bool	has_dat;			/* if true, then the data is stored in dat; otherwise: false (default) */
	public:
		/* Create a histogram based on the data x and made up of N bins; if store = true, then the data vector is stored in dat */
		Vector GetW (Vector W);
		
		Histogram( const Vector & x, const Vector & W, const int N = 20, const bool store = false );
		~Histogram(); /* Destructor */
		/* return a copy of h */
		Vector	GetHistogram( ) const ;
		void	SetHistogram( Vector & h ) const ;  /* Set? */
		/* if has_dat --> clear dat and unset has_dat */
		void	ClearData( );
		/* set new data vector and recompute the histogram */
		void	SetData( const Vector & x, const bool store = false );
		/* Return a copy of the data vector (if existing) into D and return true; otherwise: return false and delete D */
		bool	GetData( Vector &x ) const;
		/* set new number of bins and recompute; if has_dat = false --> return false (i.e. unable to change number of bins; keep old N) */
		bool	SetBinNumber( const int N );
		inline int		GetBinNumber() const { return n_bin; }
		inline bool		HasData()  const { return has_dat; }
		inline double	MinValue() const { return v_min; }
		inline double	MaxValue() const { return v_max; }
		inline double	BinWidth() const { return h_bin; }
		void		WriteDatafile( const char fn[]  ) const ;	/* Write a data file that is useable for the default GnuPlot Script (see next command) */ 
		void	WriteGnuplotScript( const char fn[] ) ;	/* Write default GnuPlot script for the creation of EPS, PNG, ... figures */
		void	WriteMatlabScript( const char fn[] ) ;	/* Write default GnuPlot script for the creation of EPS, PNG, ... figures */
};

#define _MS_DISTRIBUTION_H_
#endif /* _MS_DISTRIBUTION_H_ */

